?fastu— title: “Projeto FMS” subtitle: “AgriPredict+” author: “Octavio Deliberato Neto” date: “16 de abril de 2019” output: html_document: toc: yes editor_options: chunk_output_type: console —
O objetivo deste relatório é documentar todas as análises e etapas de modelagem que levaram ao desenvolvimento do protótipo do AgriPredict+, com dados públicos fornecidos pela Fundação MS (FMS).
Os dados brutos fornecidos pela FMS já haviam sofrido um tratamento preliminar por parte do Rodrigo Stelzer, o cientista de dados com quem se iniciou o projeto. Dele, tive acesso a um arquivo csv com dados de precipitação pluviométrica, a saber: clima_dia.csv
A primeira etapa da análise desse arquivo consiste em carregá-lo e verificar sua estrutura:
clima_dia <- read.csv("data-raw/clima_dia.csv",
stringsAsFactors = F) # arquivo do Stelzer
str(clima_dia)
'data.frame': 4018 obs. of 10 variables:
$ date : chr "2008-01-01" "2008-01-02" "2008-01-03" "2008-01-04" ...
$ AmambaÃ. : num NA NA NA NA NA NA NA NA NA NA ...
$ Anaurilandia : num 0 0 39 0 0 0 0 0 0 0 ...
$ Campo_Grande : num 0 0 3.4 2.2 4.4 28.6 0.4 29.6 0 1.2 ...
$ Dourados : num 0 0.2 1.6 28.6 12.8 0 0 0 0 0 ...
$ Ivinhema : num 8.8 44.4 0 5.4 0 0 0 0 0 0 ...
$ Maracaju : num 0.4 0 0 30.6 1.6 0.4 0.2 0 0 0 ...
$ Rio_Brilhante : num NA NA NA NA NA NA NA NA NA NA ...
$ SÃ.o_Gabriel_do_Oeste: num 19 7.6 0.4 5 0 1 0.2 30 0.8 9.4 ...
$ SidrolÃ.ndia : num NA NA NA NA NA NA NA NA NA NA ...
summary(clima_dia)
date AmambaÃ. Anaurilandia Campo_Grande
Length:4018 Min. : 0.000 Min. : 0.00 Min. : 0.000
Class :character 1st Qu.: 0.000 1st Qu.: 15.00 1st Qu.: 0.000
Mode :character Median : 0.000 Median : 22.50 Median : 0.000
Mean : 4.001 Mean : 30.51 Mean : 4.053
3rd Qu.: 1.600 3rd Qu.: 40.00 3rd Qu.: 1.800
Max. :147.600 Max. :180.00 Max. :111.000
NA's :449 NA's :3556 NA's :308
Dourados Ivinhema Maracaju Rio_Brilhante
Min. : 0.000 Min. : 0.000 Min. : 0.000 Min. : 0.00
1st Qu.: 0.000 1st Qu.: 0.000 1st Qu.: 0.000 1st Qu.: 0.00
Median : 0.000 Median : 0.000 Median : 0.000 Median : 0.00
Mean : 3.471 Mean : 3.567 Mean : 2.828 Mean : 3.99
3rd Qu.: 1.000 3rd Qu.: 0.600 3rd Qu.: 0.600 3rd Qu.: 1.20
Max. :186.000 Max. :146.000 Max. :83.400 Max. :135.00
NA's :407 NA's :372 NA's :286 NA's :369
SÃ.o_Gabriel_do_Oeste SidrolÃ.ndia
Min. : 0.000 Min. : 0.000
1st Qu.: 0.000 1st Qu.: 0.000
Median : 0.000 Median : 0.000
Mean : 1.988 Mean : 3.247
3rd Qu.: 0.400 3rd Qu.: 0.800
Max. :95.000 Max. :128.800
NA's :398 NA's :927
Nota-se o elevado número de NAs para todos os locais, notadamente Anaurilândia, e que o campo de data está representado como character, de modo que alguns ajustes se fazem então necessários:
clima_dia$date <- as.Date(clima_dia$date)
round(colSums(is.na(clima_dia))/nrow(clima_dia) * 100, 0) # % NAs
date AmambaÃ. Anaurilandia
0 11 89
Campo_Grande Dourados Ivinhema
8 10 9
Maracaju Rio_Brilhante SÃ.o_Gabriel_do_Oeste
7 9 10
SidrolÃ.ndia
23
fNames <- c("Date", "Amambai", "Anaurilandia", "Campo.Grande", "Dourados",
"Ivinhema", "Maracaju", "Rio.Brilhante", "Sao.Grabiel.do.Oeste",
"Sidrolandia")
names(clima_dia) <- fNames
clima_dia$Anaurilandia <- NULL # me livro dessa cidade, pois tem muitos NAs
clima_dia <- na.omit(clima_dia)
summary(clima_dia)
Date Amambai Campo.Grande
Min. :2008-10-01 Min. : 0.000 Min. : 0.000
1st Qu.:2010-10-30 1st Qu.: 0.000 1st Qu.: 0.000
Median :2013-06-11 Median : 0.000 Median : 0.000
Mean :2013-06-25 Mean : 4.111 Mean : 4.174
3rd Qu.:2015-10-02 3rd Qu.: 1.400 3rd Qu.: 2.400
Max. :2018-07-25 Max. :147.600 Max. :111.000
Dourados Ivinhema Maracaju Rio.Brilhante
Min. : 0.000 Min. : 0.000 Min. : 0.000 Min. : 0.000
1st Qu.: 0.000 1st Qu.: 0.000 1st Qu.: 0.000 1st Qu.: 0.000
Median : 0.000 Median : 0.000 Median : 0.000 Median : 0.000
Mean : 3.141 Mean : 3.818 Mean : 2.728 Mean : 3.992
3rd Qu.: 0.950 3rd Qu.: 0.800 3rd Qu.: 0.400 3rd Qu.: 1.200
Max. :186.000 Max. :146.000 Max. :81.400 Max. :135.000
Sao.Grabiel.do.Oeste Sidrolandia
Min. : 0.000 Min. : 0.000
1st Qu.: 0.000 1st Qu.: 0.000
Median : 0.000 Median : 0.000
Mean : 2.436 Mean : 3.363
3rd Qu.: 0.600 3rd Qu.: 0.600
Max. :95.000 Max. :128.800
A fim de agregar os dados trimestralmentre, teremos de fazer uso do pacote xts:
clima.dia.xts <- xts(clima_dia[, -1], order.by = clima_dia$Date)
clima.qrt.xts <- apply.quarterly(clima.dia.xts, mean)
ts_plot(clima.qrt.xts, type = "single")
periodicity(clima.qrt.xts)
Quarterly periodicity from 2008-12-31 to 2018-07-25
frequency(clima.qrt.xts)
[1] 1
Do gráfico anterior e de entrevistas com o Sr. Alex Melotto, presidente da FMS, parece factível agregar as precipitações dos diversos locais em uma única variável. Um mapa de correlação nos ajudará a avaliar melhor essa hipótese:
clima.cor <- cor(clima.qrt.xts)
corrplot.mixed(clima.cor)
ts_plot(clima.qrt.xts[, c("Rio.Brilhante", "Ivinhema")])
À exceção de São Gabriel do Oeste, fica evidente a lucidez dessa proposta, i.e. agregar as precipitações dos locais, que é o que se faz a seguir:
clima.qrt <- as.data.frame(clima.qrt.xts)
clima.qrt$Date <- as.Date(rownames(clima.qrt))
rownames(clima.qrt) <- NULL
clima.qrt$quarter <- quarter(clima.qrt$Date)
clima.qrt$harvest <- year(clima.qrt$Date)
clima.qrt <- clima.qrt %>% select(harvest, quarter, everything(), -Date)
clima.qrt$rain.mm <- rowMeans(clima.qrt[, 3:ncol(clima.qrt)])
clima.qrt <- clima.qrt %>% select(harvest, quarter, rain.mm)
str(clima.qrt)
'data.frame': 36 obs. of 3 variables:
$ harvest: num 2008 2009 2009 2009 2009 ...
$ quarter: int 4 1 2 3 4 1 2 3 4 4 ...
$ rain.mm: num 3.55 5.23 1.49 2.94 5.2 ...
rio::export(x = clima.qrt, file = "clima_qrt.csv")
Finalmente, por solicitação do próprio Sr. Alex Melotto, seria melhor que os dados trimestrais de precipitação fossem traduzidos em classes, e.g. “sem chuva”, “pouca chuva”, “chuva média”, “muita chuva”, em vez de números.
# Níveis de precipitação
rain.level <- cut2(clima.qrt$rain.mm, g = 4)
levels(rain.level) <- c("norain", "low", "med", "high")
clima.qrt$rain.level <- as.character(rain.level)
clima.qrt$rain <- paste0('q', clima.qrt$quarter, '.', clima.qrt$rain.level)
clima.qrt$rain <- factor(clima.qrt$rain, levels = c(
'q1.norain', 'q1.low', 'q1.med', 'q1.high',
'q2.norain', 'q2.low', 'q2.med', 'q2.high',
'q3.norain', 'q3.low', 'q3.med', 'q3.high',
'q4.norain', 'q4.low', 'q4.med', 'q4.high'))
# Um pouco de feature engineering...
clima.qrt$rain.level <- NULL
dummies <- dummy_cols(clima.qrt)
dummies <- dummies %>% select(harvest,
rain_q1.low, rain_q1.med, rain_q1.high,
rain_q2.norain, rain_q2.low, rain_q2.med,
rain_q3.norain, rain_q3.low,
rain_q4.low, rain_q4.med, rain_q4.high)
clima.qrt.final <- dummies %>%
group_by(harvest) %>% summarise_if(is.numeric, max)
str(clima.qrt.final)
Classes 'tbl_df', 'tbl' and 'data.frame': 11 obs. of 12 variables:
$ harvest : num 2008 2009 2010 2011 2012 ...
$ rain_q1.low : num 0 0 0 0 0 0 0 1 0 0 ...
$ rain_q1.med : num 0 0 0 0 1 0 0 0 1 1 ...
$ rain_q1.high : num 0 1 1 0 0 1 1 0 0 0 ...
$ rain_q2.norain: num 0 1 1 0 0 0 0 0 1 0 ...
$ rain_q2.low : num 0 0 0 0 0 1 1 0 0 1 ...
$ rain_q2.med : num 0 0 0 0 1 0 0 0 0 0 ...
$ rain_q3.norain: num 0 0 1 0 1 1 0 0 0 1 ...
$ rain_q3.low : num 0 1 0 0 0 0 1 1 1 0 ...
$ rain_q4.low : num 0 0 1 0 0 0 0 0 0 0 ...
$ rain_q4.med : num 1 0 0 0 0 1 1 1 1 0 ...
$ rain_q4.high : num 0 1 0 1 1 0 0 0 0 1 ...
Finalmente, salvo o arquivo resultante:
write_csv(clima.qrt.final, "clima.csv")
Os dados brutos de produtividade das culturas de soja já haviam sofrido um tratamento preliminar por parte do Rodrigo Stelzer, o cientista de dados com quem se iniciou o projeto. Dele, tive acesso a um arquivo csv com esses dados, a saber: crop.csv
A primeira etapa da análise desse arquivo consiste em carregá-lo e verificar sua estrutura:
crop <- read.csv("data-raw/crop.csv", encoding="UTF-8",
stringsAsFactors = F) # arquivo do Stelzer
str(crop)
'data.frame': 14436 obs. of 6 variables:
$ id : int 0 1 2 3 4 5 6 7 8 9 ...
$ material : chr "fts campo mourão rr" "na 5909 rg" "brs 282" "brs 284" ...
$ harvest : int 2008 2008 2008 2008 2008 2008 2008 2008 2008 2008 ...
$ site : chr "antônio joão" "antônio joão" "antônio joão" "antônio joão" ...
$ productivity: num 66 65 63.6 63.5 59.9 58.4 57.8 57.3 56.6 56.2 ...
$ species : chr "soy" "soy" "soy" "soy" ...
summary(crop)
id material harvest site
Min. : 0 Length:14436 Min. :2008 Length:14436
1st Qu.:1804 Class :character 1st Qu.:2010 Class :character
Median :3608 Mode :character Median :2013 Mode :character
Mean :3612 Mean :2013
3rd Qu.:5413 3rd Qu.:2015
Max. :7444 Max. :2018
productivity species
Min. : 2.60 Length:14436
1st Qu.: 58.60 Class :character
Median : 72.80 Mode :character
Mean : 84.41
3rd Qu.:114.10
Max. :198.60
length(unique(crop$material))
[1] 835
barplot(table(crop$material))
Da análise anterior, nota-se um grande número de sementes distintas, cujas classes estão nitidamente desbalanceadas, o que pode ser um problema para o treinamento dos algoritmos de machine learning.
A seguir, mais algumas análises e gráficos de exploração:
crop$harvest <- factor(crop$harvest)
crop$id <- NULL
boxplot(productivity ~ harvest, data = subset(crop, species == 'soy'),
main = 'Produtividade da Soja')
boxplot(productivity ~ harvest, data = subset(crop, species == 'corn'),
main = 'Produtividade do Milho')
table(crop$harvest)/nrow(crop)*100
2008 2009 2010 2011 2012 2013 2014
1.780272 11.277362 12.835966 10.092823 7.460515 10.944860 11.360488
2015 2016 2017 2018
9.843447 12.524245 9.434746 2.445276
cumsum(table(crop$harvest)/nrow(crop)*100) # até 2015 para treinamento, resto para teste
2008 2009 2010 2011 2012 2013
1.780272 13.057634 25.893599 35.986423 43.446938 54.391798
2014 2015 2016 2017 2018
65.752286 75.595733 88.119978 97.554724 100.000000
g <- ggplot(data = subset(crop, species == 'soy'),
aes(x = harvest, y = productivity)) + geom_boxplot() +
facet_wrap(~ site, ncol = 4) + labs(title = 'Produtividade da Soja')
g + theme(axis.text.x = element_text(angle = 90, hjust = 1))
g1 <- ggplot(data = subset(crop, species == 'corn'),
aes(x = harvest, y = productivity)) + geom_boxplot() +
facet_wrap(~ site, ncol = 4) + labs(title = 'Produtividade do Milho')
g1 + theme(axis.text.x = element_text(angle = 90, hjust = 1))
crop.sum <- crop %>% group_by(species, site) %>%
summarise(avgProd = mean(productivity), sdProd = sd(productivity)) %>%
arrange(desc(avgProd))
crop.sum
# A tibble: 32 x 4
# Groups: species [2]
species site avgProd sdProd
<chr> <chr> <dbl> <dbl>
1 corn sidrolândia 128. 26.6
2 corn bonito 126. 16.7
3 corn amambai 118. 26.2
4 corn são gabriel do oeste 116. 22.7
5 corn rio brilhante 115. 24.8
6 corn maracaju 113. 22.1
7 corn antônio joão 110. 10.7
8 corn dourados 108. 38.8
9 corn naviraí 108. 20.7
10 corn deodápolis 99.9 23.4
# ... with 22 more rows
g2 <- ggplot(crop, aes(x = productivity, fill = species)) +
geom_density(alpha = 0.5)
g2
summary(crop.sum$avgProd)
Min. 1st Qu. Median Mean 3rd Qu. Max.
34.07 55.94 61.65 75.30 107.88 128.45
g3 <- ggplot(data = crop.sum, aes(x = site, y = avgProd, color = species)) +
geom_point(size = 2)
g3 + theme(axis.text.x = element_text(angle = 90, hjust = 1))
Dos gráficos, nota-se que há variação significativa entre os lugares, o ano de plantio e o tipo de cultura. De alguma maneira, todas essas variáveis deverão, então, constar do modelo final.
A seguir, unem-se os dados climáticos com os de produtividade, a fim de formar o dataset que será usado na etapa de modelagem:
clima <- read_csv("clima.csv")
crop.df <- merge(crop, clima, by = 'harvest')
crop.df$material <- str_replace_all(crop.df$material, "[^[:alnum:]]", "")
seeds <- substr(crop.df$material, 1, 2)
tbl.seeds <- table(seeds)
sum(tbl.seeds/sum(tbl.seeds)*100) # check sum=100
[1] 100
barplot(tbl.seeds/sum(tbl.seeds)*100) # types of materials
crop.df$site <- str_replace_all(crop.df$site, "[^[:alnum:]]", "")
crop.df$site <- str_replace_all(crop.df$site, "[[:punct:]]", "")
crop.df$material <- seeds
crop.df <- crop.df %>% select(material, site, species, contains('rain_q'),
productivity)
write_csv(crop.df, 'crop_df.csv')
O arquivo final, de nome crop_df.csv, contém os dados necessários à etapa de modelagem matemática da produtividade.
DT::datatable(crop.df)
library(caret)
library(ranger)
crop_df <- read_csv("crop_df.csv")
crop_df$material <- substr(crop_df$material, 1, 1) # try to balance things out
barplot(table(crop_df$material)) # still not balanced, but let's give it a shot
crop_df[, 1:14] <- lapply(crop_df[, 1:14], as.factor)
## Not run (takes a lot of time):
## Random Forest
# ctrl <- trainControl(method = "cv", selectionFunction = 'oneSE')
# grid <- expand.grid(.mtry = c(6, 8, 10), .splitrule = c("extratrees"),
# .min.node.size = c(1, 3, 5))
# set.seed(300)
# dim(crop_df)
# m.rf <- train(productivity ~ ., data = crop_df, method = "ranger",
# metric = "RMSE", trControl = ctrl, tuneGrid = grid)
# m.rf
# yhat.rf <- predict(m.rf, crop_df)
# plot(crop_df$productivity, yhat.rf)
# abline(lm(yhat.rf ~ crop_df$productivity), col = 'darkred', lwd = 2)
# cor(crop_df$productivity, yhat.rf)^2
# saveRDS(m.rf, 'm_rf.rds')
## End (Not run)
# In order to save the time required to train the RF model
m.rf <- readRDS('m_rf.rds')
m.rf
Random Forest
14436 samples
14 predictor
No pre-processing
Resampling: Cross-Validated (10 fold)
Summary of sample sizes: 12993, 12993, 12992, 12993, 12991, 12994, ...
Resampling results across tuning parameters:
mtry min.node.size RMSE Rsquared MAE
6 1 14.71343 0.8417227 10.933627
6 3 14.66929 0.8426394 10.907607
6 5 14.63400 0.8432560 10.872263
8 1 13.18050 0.8636927 9.683304
8 3 13.18721 0.8637343 9.695281
8 5 13.20974 0.8634428 9.710814
10 1 12.51664 0.8729141 9.101445
10 3 12.50263 0.8731346 9.094314
10 5 12.54313 0.8725464 9.126912
Tuning parameter 'splitrule' was held constant at a value of extratrees
RMSE was used to select the optimal model using the one SE rule.
The final values used for the model were mtry = 10, splitrule =
extratrees and min.node.size = 1.
yhat.rf <- predict(m.rf, crop_df)
plot(crop_df$productivity, yhat.rf,
main = "Random Forest",
xlab = "Produtividade, sacas / ha",
ylab = "Produtividade calculada, sacas / ha")
abline(lm(yhat.rf ~ crop_df$productivity), col = 'darkred', lwd = 2)
r2 <- cor(crop_df$productivity, yhat.rf)^2
text(25, 140, paste0("R2 = ", round(r2, 3)))